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We  present  an  asymptotic  study  on  the  phase  diagram  of  two-dimensional  nematic  liquid  crystal 
polymer  monolayers  with  the  Onsager  intermolcular  potential.  In  contrast  to  the  case  of  Maier- 
Saupe  interaction  potential  where  there  is  only  one  nematic  branch,  our  analysis  reveals  that  there 
are  infinite  many  nematic  branches  in  the  case  of  the  Onsager  interaction  potential.  An  asymptotic 
expression  is  derived  for  each  nematic  branch.  For  small  polymer  concentration  the  isotropic  branch 
is  the  only  equilibrium  state.  As  the  polymer  concentration  is  increased,  nematic  branches  appear 
one  by  one,  starting  with  the  first  nematic  branch.  The  polymer  orientation  distribution  of  the  first 
nematic  branch  has  a  two  fold  rotational  symmetry,  the  second  branch  has  a  four  fold  rotational 
symmetry,  the  third  branch  has  a  six  fold  rotational  symmetry,  and  so  on.  To  determine  the  stability 
of  these  nematic  branches,  we  derive  an  asymptotic  expression  of  free  energy  for  each  nematic 
branch.  We  find  that  free  energies  of  all  nematic  branches  are  lower  than  that  of  the  isotropic 
state,  and  the  first  nematic  branch  has  the  lowest  free  energy  among  all  branches.  To  further 
investigate  the  stability  and  meta-stability,  we  carry  out  asymptotic  analysis  of  the  free  energy  when 
each  nematic  state  is  perturbed.  We  conclude  that  (1)  the  isotropic  branch  is  stable  until  the  first 
nematic  branch  appears,  (2)  the  first  nematic  branch  is  stable,  and  (3)  the  isotropic  branch  (after 
the  appearance  of  the  first  nematic  branch)  and  all  other  nematic  branches  are  unstable  when 
perturbed  by  the  leading  Fourier  mode  in  the  first  nematic  branch.  We  also  present  a  spectrum 
numerical  method  for  calculating  nematic  branches  and  free  energies.  The  spectrum  method  yields 
results  that  are  accurate  up  to  the  computer  precision.  All  of  asymptotic  results  are  confirmed  by 
numerical  results  obtained  with  the  spectrum  method. 
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1.  INTRODUCTION 

Microstructures  play  a  significant  role  in  understanding 
the  phase  transition  for  dispersions  of  rodlike  macro¬ 
molecules  in  liquid  crystal  polymers  with  or  without  exter¬ 
nal  field.  In  1949  Onsager  developed  the  first  theory  that 
demonstrated  an  isotropic-nematic  (IN)  phase  transition 
due  to  steric  interactions  alone.33  His  theory  was  based 
on  a  virial  expansion  which  ended  up  with  a  compli¬ 
cated  nonlinear  integral  equation  for  the  orientation  distri¬ 
bution  function.  By  using  a  nematic  potential  (“Onsager 
potential’')  and  choosing  an  appropriate  trial  function  for 
the  orientation  distribution  function,  Onsager  was  able 
to  argue  that  when  the  polymer  concentration  is  high 
enough,  there  is  a  transition  from  a  uniform  isotropic 
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state  to  an  orientationally  ordered  prolate  nematic  state. 
The  Onsager  theory  is  valid  for  dilute  polymer  solutions 
because  the  virial  expansion  is  terminated  at  the  second 
virial  coefficients. 

Many  theories  have  been  formulated  after  Onsager’s 
work.  For  an  introduction  to  the  theoretical  advances  of  liq¬ 
uid  crystalline  polymers  we  refer  the  reader  to  the  review 
paper  by  Rey  and  Denn.34  Among  many  theories  the  Doi- 
Hess  kinetic  model6, 20  is  gaining  popularity  in  recent  years 
for  describing  various  properties  of  liquid  crystalline  poly¬ 
mers  (LCPs)  in  a  solvent.  This  well-known  kinetic  the¬ 
ory  was  first  developed  by  Doi  and  Edwards  to  describe 
spatially  homogeneous  flows  of  rodlike  LCPs.  The  sig¬ 
nificance  of  the  Doi-Hess  theory  lies  on  the  fact  that 
it  models  the  ensemble  of  rodlike  macromolecules  using 
a  meso-scale  probability  density  function.  Thus,  it  can 
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accommodate  the  isotropic  state,  the  nematic  state,  and 
the  phase  transition  between  them.  In  the  Doi-Hess  kinetic 
model,  the  intermolecular  interaction  is  usually  modeled 
by  the  Maier-Saupe  potential  because  of  the  simplicity  of 
the  Maier-Saupe  potential.  The  Maier-Saupe  potential  is  an 
approximation  to  the  Onsager  potential.  In  the  Maier-Saupe 
potential,  the  effect  of  the  polymer  ensemble  on  one  poly¬ 
mer  rod  is  represented  by  the  second  moment  of  the  prob¬ 
ability  density.  Thus,  with  the  Maier-Saupe  potential,  an 
equilibrium  state  is  completely  determined  by  the  second 
moment,  which  motivates  the  idea  of  following  the  evolu¬ 
tion  of  only  the  second  moment  with  various  closures.  A  lot 
of  studies  have  been  devoted  to  the  Doi-type  kinetic  theory 
and  its  approximation  with  various  closures  based  on  the 
Maier-Saupe  potential  and  the  Onsager  potential  (for  exam¬ 
ple,  Refs.  [3-5,  8-17,  19,  23-28,  30,  35-38,  41,  43-49]). 
However,  almost  nothing  can  be  found  in  the  literature  for 
the  theoretical  analysis  of  the  Doi-Hess  kinetic  model  with 
the  Onsager  interaction  potential. 

Recently  the  authors  carried  out  a  numerical  study  on 
the  equilibrium  states  of  rigid  polymer  rods  with  the 
Onsager  potential.40  The  primary  purpose  of  the  current 
paper  is  to  extend  the  work  in  Ref.  [40]  and  conduct 
a  theoretical  study  on  the  full  Doi-Hess  kinetic  theory 
with  the  Onsager  interaction  potential.  We  confine  our 
attention  to  the  two-dimensional  liquid  crystal  polymers 
which  is  physically  motivated  by  monolayer  films  in 
microelectronics  and  optoelectronics.  Rigid,  rodlike  liq¬ 
uid  crystalline  polymers  have  been  of  considerable  interest 
in  thin  films  because  of  their  stability,  nonlinear  optical 
properties,  and  role  as  alignment  layers  for  liquid  crys¬ 
tal  displays.29  Experiments  have  shown  that  the  dynamic 
behavior  of  monolayer  nematic  polymers  can  be  quanti¬ 
tatively  captured  by  the  simple  molecular  model  of  Doi 
and  Hess  for  both  extensional  and  simple  shear  flows.29' 32 
Numerical  and  theoretical  studies  on  the  two-dimensional 
liquid  crystal  polymers  have  also  been  conducted 
by  many  researchers. ',2' 7' 18'21_23,26_28'30'31'39'42,44' 47’ 49  Our 
work  extends  these  previous  papers  to  the  Doi-Hess  theory 
with  the  Onsager  intermolecular  potential. 

This  paper  proceeds  as  follows.  First,  the  mathematical 
framework  of  the  two-dimensional  Doi-Hess  model  with 
the  Onsager  potential  is  described  in  Section  2,  and  the  con¬ 
nection  between  the  Onsager  potential  and  the  Maier-Saupe 
potential  is  briefly  discussed.  In  Section  3,  an  asymptotic 
analysis  is  carried  out  to  calculate  the  nematic  branches 
of  polymers  with  the  Onsager  potential.  It  turns  out  that 
there  are  infinitely  many  nematic  branches  for  polymers 
with  the  Onsager  potential.  A  spectrum  method  is  described 
in  Section  4.  The  spectrum  method  is  capable  of  reduc¬ 
ing  numerical  errors  to  the  computer  precision.  Numerical 
results  obtained  with  the  spectrum  method  are  treated  as 
the  exact  solution  and  are  used  to  confirm  the  asymptotic 
results.  In  Section  5,  the  stability  of  the  nematic  branches 
is  analyzed  using  free  energy  calculations.  Finally,  in 
Section  6,  we  summarize  the  results  obtained  in  this  paper. 


2.  MATHEMATICAL  FORMULATION  OF  THE 
TWO  DIMENSIONAL  DOI-HESS  KINETIC 
MODEL  WITH  THE  ONSAGER  POTENTIAL 


We  briefly  review  the  two-dimensional  mathematical  for¬ 
mulation  of  the  Doi-Hess  kinetic  theory  for  rigid  rodlike 
nematogenic  molecules  immersed  in  a  viscous  solvent.6, 20 
In  this  theory,  one  considers  a  test  molecule  in  a  sea 
of  others  and  models  the  interactions  between  molecules 
by  a  mean-field  potential.  The  ensemble  of  polymer  rods 
is  described  by  an  orientation  probability  density  func¬ 
tion  (PDF).  In  the  absence  of  flow,  the  orientation  prob¬ 
ability  density  function  is  governed  by  the  Smoluchowski 
equation: 

dp  d  (  1  dV  dp\ 

—  =D — •( - P+—  (1) 

dt  du  \kBT  du  duj 

where  p(u,  t)  is  the  probability  density  function  that  a  poly¬ 
mer  rod  has  orientation  u  at  time  t,  u  is  a  unit  vector 
denoting  the  polymer  rod  orientation,  d/du  is  the  gradient 
operator  on  the  unit  sphere,6  D  the  rotational  diffusivity 
constant,  kB  the  Boltzmann  constant,  T  the  absolute  tem¬ 
perature  and  V  the  mean-field  potential  representing  the 
interactions  between  polymer  rods. 

In  the  two-dimensional  space,  the  orientation  of  a  poly¬ 
mer  rod  can  be  described  by  an  angle  0.  The  orientational 
vector  of  a  polymer  rod  is  given  by  u  =  (cos  0 ,  sin  9).  The 
Onsager  intermolecular  potential  describes  the  interaction 
between  the  test  polymer  rod  and  all  other  polymer  rods. 
The  Onsager  potential  has  the  form 

U(9)  ee  ^nTr  =  b  r  | sin (0  9)\p(0)  dO  (2) 

kb  I  J  o 

where  b  is  the  strength  of  the  Onsager  potential  which 
is  proportional  to  the  normalized  polymer  concentration 
and  is  inversely  proportional  to  the  absolute  temperature, 
and  p(0)  the  probability  density  function  of  the  orientation 
angle  6.  For  the  convenience  of  discussions  below,  we  use 
U{9)  to  denote  the  Onsager  potential  normalized  by  the 
factor  kBT.  It  should  be  pointed  out  that  the  Maier-Saupe 
intermolecular  potential  has  a  similar  form  as  the  Onsager 
potential  with  the  kernel  |sin(0  —  0)\  in  (2)  being  replaced 
by  |sin(0—  9)\2  =  (1  —  cos 2(0  —  9))/2.  The  Maier-Saupe 
potential  has  the  form 


t/(0)  .  hf  =  i,  f 

kvl  J  o 


^  1-  cos  2(6-0) 


kBT 
b  b 
2  ~~  2 


p(9)  dO 


n2l T  ~  ~  ~ 

cos  2 9  /  cos  2 6p(6)d0 

n2.1T  ~  ~  J 

+  sin20/  sin2  9p(0)d9 
do 


(3) 


In  the  Maier-Saupe  potential  (3),  the  effect  of  p{9)  is  com¬ 
pletely  specified  by  two  quantities  cos  29  p(9)d  9  and 
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f'77  sin29p(9)d0.  In  other  words,  if  these  two  quantities 
are  given,  then  the  Maier-Saupe  potential  is  completely 
determined  without  knowing  the  probability  density  p{9). 
In  contrast,  the  Onsager  potential  (2)  depends  on  the  whole 
function  of  the  probability  density  p(9).  This  property  of 
the  Onsager  potential  causes  significant  mathematical  dif¬ 
ficulties  for  theoretical  analysis. 

In  the  current  study  of  the  two  dimensional  Doi-Hess 
kinetic  model,  we  take  the  advantage  that  the  probability 
density  p(6)  is  a  periodic  function.  We  consider  the  Fourier 
series  expansion  of  the  probability  density  function  p(6): 

1 

p(8)  =  - - \-y2(akcoslc0  +  bksinkd)  (4) 

2  77 


where  the  constant  term  being  1/(277)  is  dictated  by  the 
fact  that  the  total  probability  is  one:  f~n  p (9) d 9  —  1. 
Below,  we  will  use  the  Fourier  series  of  p(9)  to  calculate 
the  Fourier  expansion  of  the  normalized  Onsager  potential 
U{0).  We  start  with  some  mathematical  preparations: 

o2tt 

/  \sin9\cosk9d9 
J o 


sindcos  kOdO  — 
sinffcos  kOdO  — 

IT 

sindcos  kOdO 


n2lT 

/  sindcos  kddO 

J  7T 

f  sin(2ir  —  0)cos(2kTT  —  k9)d8 
Jo 


(  2  2 
J  k+  1  k- 1 


4 

- ,  for  k  even 

k2-\ 


0, 


for  k  odd 


(5) 


/» 2  77"  /»  77"  /» 2  77" 

/  |sin0|sink0(f0  =  /  sinflsin kddO—  /  sinOsinkddO 

j  0  J  0  "  TT 

=  /  sinOsinkOdO 
A) 

—  [  sin(2TT  —  0)sin(2kir  —  k0)d0  =  0 
Jo 

(6) 

With  the  help  of  (4),  (5)  and  (6),  we  write  the  normalized 
Onsager  potential  as 


=  b 

o2tt 

A) 

|sin(0 

-9)\p(8)d0 

n2,TT- 

-0 

=  b 

/ 

sin  9\p(9  +  9)d9 

'0-0 

nllT 

=  b 

L 

sin  9\p(9  +  9)d9 

T 

"  i 

=  b 

l 

|sin0| 

—  +  Y,(akcosk(Q+0) 
L27T  k= 1 

+  bksink(9  +  9)) 
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=  b 


|sin0| 


1 

277 


oo 

^2(ak  cos  k9 cos  k9 

k=  1 


—  ak  sinkd  sin  k8) 


d9 
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+b  |sin0|  y ~^{bk  sinkiicoskfl 
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+bk  cos  k9  sin  k9)d9 
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= - 4  i>-V 

277  ~l  L  (2j)2  —  1 
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cos  2 jd 
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(2/)2  —  1 

For  mathematical  convenience,  we  introduce 

4  bb0l 


4b  a2j 


d2 ;  — 


V 


(2 7)2-r  22  (2;)2-i 


(7) 

(8) 


We  use  c2j  and  d2j  as  unknowns,  and  treat  ak  and  bk  as 
functions  of  c2;  and  d2)  ■  Since  the  intermolecular  interac¬ 
tion  potential  drives  polymer  rods  only  through  its  deriva¬ 
tive,  a  constant  term  in  the  potential  will  not  change  the 
dynamics.  We  drop  the  constant  term  in  (7)  and  use  the 
new  variables  c2j  and  d2j  to  write  the  Onsager  potential  as 


U(8)  =  -  XXC2 j cos  2J0  +  dij sin  2J°)  (9) 

j=  i 

The  equilibrium  probability  density  of  Smoluchowski 
equation  (1)  is  given  by  the  Boltzmann  distribution: 

P (0)  =  ^exp  [-1/(0)] 

=  \  exP  Y,(c2j cos  2je  +  d2j  sin  2j9)^j  (10) 
where  Z  is  the  partition  function  given  by 


Z  = 


=  f  exp[—  U(8)]d0 
Jo 

=  J  exp  ^ XXc2j cos 2j8  +  d2j  sin  2 j8)'j d9  (11) 


Recall  that  we  started  with  the  Fourier  series  expansion  of 
p(8)  with  coefficients  a,  and  I? The  Fourier  coefficients 
with  even  indices,  a2k  and  b2k,  can  be  expressed  in  terms 
of  c2k  and  d2k  from  (8). 

a2k  =  2m)2-l]c2k,  b2k  —  2\(2k)2  —  V)d2k  (12) 

Below  we  use  (10)  to  show  that  the  Fourier  coefficients 
of  p{9)  with  odd  indices  are  zero.  The  probability  density 
p(9)  given  in  (10)  satisfies 


p(0  +  77)  =  p(0)  (13) 
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With  the  help  of  (13),  we  have 


*2k+l 


1  r 2'n 

—  cos[(2A'+  l)9]p(9)d9 

TT  Jo 

—  [  cos[(2k+ l)9]p(9)d9 
TT  Jo 

1  r2n 

H - /  cos[(2k+l)9]p(9)d9 

TT  J  i t 

—  [  cos[(2&+  l)6]p(9)dd 
TT  Jo 

—  f  cos[(2fc  +  1)0+  7r]p(0  +  Tr)dd  =  0 
tt  Jo 


+  — 
77 


(14) 


Similarly,  one  can  verify  that 
1  r2 * 

b2k+1  =  -  /  sin[(2 k+  1  )9]p(0)d9  =  0  (15) 

TT  Jo 

Next  we  show  that  U (9)  can  be  made  an  even  function  of 
9  by  shifting  the  coordinate  system.  That  is,  we  can  make 
b2k  —  0  for  all  A:  by  a  single  shifting.  We  start  by  using 
relation  (12)  to  express  function  p(9)  in  terms  of  function 
U(9)  and  its  derivatives. 

P(.0)=^-  +  ^r[U”(.8)  +  um  (16) 

277  4 b 

The  equilibrium  probability  density  function  p{9)  is  also 
linked  to  the  normalized  Onsager  potential  U  (6)  through 
the  Boltzmann  relation 

p(9)=  ^exp[-U(d)\  (17) 

Combining  (16)  and  (17),  we  obtain  a  differential  equation 
for  U(9). 

C"(0)  =  -£/(0)+^exp[-f/(0)]-^  (18) 

Z  277 

(18)  is  a  second  order  differential  equation.  To  determine 
the  solution  of  a  second  order  differential  equation,  we  need 
to  specify  two  initial  conditions.  Let  90  be  the  location  at 
which  p{9)  attains  its  global  maximum.  It  follows  immedi¬ 
ately  from  (17)  that  the  normalized  Onsager  potential  U(9) 
attains  its  global  minimum  at  90.  Without  loss  of  general¬ 
ity,  we  assume  90  —  0  (one  can  always  shift  the  coordinate 
system  to  make  0O  —  0).  At  the  location  of  the  global  max¬ 
imum,  we  have 

U  (0)  =  UQ,  U\ 0)  =  0  (19) 

We  use  (19)  as  the  initial  conditions  for  differential  equa¬ 
tion  (18).  It  is  straightforward  to  verify  that  both  differential 
equation  (18)  and  initial  conditions  (19)  are  invariant  under 
the  change  of  variable  9-^—9.  Therefore,  we  conclude 
that  U (9)  is  an  even  function 

U(9)  =  U(—9)  (20) 


Consequently,  from  (9)  we  have  d2k  =  0  and  U  ( 0 )  is  sim¬ 
plified  to 

00 

U(9)  =  ~Y,c2kCos2k9  (21) 

k=  1 

The  next  task  is  to  derive  a  (nonlinear)  system  for  Fourier 
coefficients  c2k.  Relation  (12)  provides  one  set  of  equations 
for  Fourier  coefficients  a2k  and  c2k.  To  derive  another  set 
of  equations  for  a2k  and  c2k,  we  calculate  the  Fourier  coef¬ 
ficients  of  p{9)  using  the  expression  of  p(9)  given  in  (10) 


1  r2n 

a2k  =  —  \  cos  2  k9p{9)d9 
TT  Jo 

1  cos2A'0exp[—  U(9)]d9 

~  77  Cexp[-U(9)\d9 

_  J_  IT  cos  2k 9  exp  K”  1  c2j  cos  2 id]  d9 
77  /(T  eXP  1  c2i  cos  2W] d0 

Combining  (12)  and  (22),  we  arrive  at  a  nonlinear  system 
for  coefficients  { c2k } . 


/02lT  cos  2k9  exp  [  j  c2i  cos  2i9]  d9 

$T  exP  [T,T=  l  c2i cos  2'0]  db 


TT 

4 b 


[(2A')2-1]c2A. 


(23) 


Since  the  Boltzmann  relation  establishes  a  one-to-one  cor¬ 
respondence  between  the  probability  density  p(9)  and  the 
Onsager  potential  U(9)  (if  we  set  the  constant  term  to 
zero  in  the  Fourier  expansion  of  the  Onsager  potential),  we 
only  need  to  solve  for  the  Onsager  potential.  In  a  prop¬ 
erly  selected  coordinate  system,  the  Onsager  potential  is  an 
even  function  and  its  Fourier  coefficients  satisfy  c2k+1  =  0 
and  satisfy  nonlinear  system  (23).  The  nonlinear  system 
(23)  is  the  key  governing  equation  in  our  study.  Below  we 
will  solve  it  first  asymptotically  and  then  numerically  with 
a  very  accurate  numerical  method. 


3.  PHASE  DIAGRAM  OF  NEMATIC 
POLYMERS  WITH  THE  ONSAGER 
POTENTIAL 

We  now  solve  nonlinear  system  (23)  to  obtain  the  phase 
diagram  of  nematic  polymers  with  the  Onsager  interaction 
potential.  We  carry  out  asymptotic  analysis  for  the  case 
where  the  amplitude  of  the  normalized  Onsager  potential  is 
small.  Since  the  probability  density  and  the  Onsager  poten¬ 
tial  are  related  by  the  Boltzmann  distribution,  small  ampli¬ 
tude  in  the  Onsager  potential  implies  that  the  probability 
density  is  not  far  from  being  a  constant.  That  is,  we  do 
asymptotic  analysis  on  the  nematic  states  that  are  close  to 
the  isotropic  state.  Mathematically,  we  assume  the  Fourier 
coefficients  of  the  Onsager  potential  have  the  expansions 

c2j  =  ea2j  +  s2/ 32j  +  0(s3)  (24) 

where  £  is  a  small  parameter,  and  «2/  and  p2j  are  coef¬ 
ficients  independent  of  e.  At  this  point,  the  meaning  of  e 
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is  not  precisely  defined  yet.  Later  we  will  re-define  e  pre¬ 
cisely.  Expanding  the  numerator  and  denominator  on  the 
left  side  of  (23)  in  powers  of  b  leads  to 


nllT 

00 

/  cos 2k 0  exp 
Jo 

J2(ea2i+e2f32i)cos2i8 

_i=  i 

n2.Tr 

=  /  cos2 kO 

J  o 


dd 


1  +  J2(sa2i  +  e2/32 /)  cos2  id 
1=1 

+  ^e2  (j2a2icos2ie\  +  0(e3) 


dd 


=  (sa2  k  +  E2p2k)ir 


J  2tt  /  00  \ 

+  -e2J  cos2kdy22a2icos2i0J  dd+0(s 3) 


(25) 


and 


n2ir 

00 

/  exp 

XXea2i  +  fi2/^2i)  COS  2id 

Jo 

_i= 1 

dO 


1  +  'y'](sa0,  +  e2/32i)  cos  2  id  +  0(e2) 

i=t 


dO 


—  2tt[1  +  0(e2)] 


(26) 


respectively.  In  general,  we  expect  the  Onsager  potential 
U{6)  changes  with  the  normalized  polymer  concentra¬ 
tion  /;,  which  suggests  that  b  is  not  independent  of  e.  We 
assume  b  has  the  expansion 

b  =  bm[l+sb(1)  +  0(s2)]  (27) 

Substituting  these  results  into  Eq.  (23),  and  first  looking 
at  only  the  coefficients  of  the  first  order  (0(e))  terms,  we 
obtain 

£^  =  £4 ^[(2A')2-1]“-  (») 

which  in  turn  yields 

alk{l-^[(2k)2-l]}=0  (29) 

Below  we  consider  possible  solutions  of  (29)  in  several 
cases. 

•  Case  I :  If  a2  ^  0,  then  in  (29)  for  k  =  1  we  must  have 
1-2^[22-1]  =  0 

which  yields  =  ((22  —  l)/2)t7  =  (3/2)tt.  Conse¬ 
quently,  for  k  >  1,  we  have  1  —  (7r/2f>(0))[(2A:)2  —  1]  = 
(22  -  (2k)z)/(22  -  1)  <  0.  It  follows  from  (29)  that  a2k  =  0 
for  k  >  1 .  Recall  that  the  meaning  of  e  in  the  expansion  is 
not  precisely  defined  yet.  When  a2  /  0,  we  re-define  c2  = 
a2s  +  /32£2  +  0(s3)  as  the  new  s,  which  is  equivalent  to 


setting  a2  =  1  and  f32  =  0.  Substituting  a2=  1  and  a2k  =  0 
for  k  >  1  into  (25),  we  have 


RHS  of  (25) 


=  (ea2k  +  E2l32k)'TT  + 


k  =  2 

+  0(e 3)  (30) 

otherwise 


In  the  calculation  of  coefficients  a2k  above,  we  keep  only 
the  first  order  terms  in  Eq.  (23).  Now  in  the  calculation  of 
coefficients  f3lk,  we  keep  both  the  first  order  and  the  second 
order  terms.  We  first  calculate  /;(l).  For  k  =  1,  Eq.  (23) 
becomes 


e  +  0(e3)  e 

1  +  0(£2)  “  l  +  eM9  +  0(e2) 

which  yields  =  0.  For  k  =  2,  (23)  provides  us  an  equa¬ 
tion  for  (3t: 

e2/34  +  s21/4+0(b3)  _  42-l  e2/34 

1  +0(s2)  _  22  -  l'  1  +  0(s2) 

which  leads  to  /34  =  1/16.  For  k  >  2,  Eq.  (23)  takes  the 
form 


£2j32A.  +  Q(g3)_  (2A-)2-l  E2(32k 

1  +  0(s2)  3  1  +  0(s2) 

which  gives  us  f32k  =  0  for  k  >  2.  Therefore,  the  nematic 
state  that  has  a  small  non-zero  coefficient  for  cos  26  in  the 
Fourier  expansion  of  U (9)  is  described  by 

E2 

U  (9)  =  — scos20 - cos40+  0(e3) 

16 

b=^TT(l  +  0(E2))  (31) 

To  calculate  the  coefficient  of  the  second  order  term  in  the 
expansion  of  b,  we  keep  one  more  term  in  the  expansions 
of  (25)  and  (26): 


5  77 

RHS  of  (25)=  stt  +  s3—  +  0(e4),  k  =  1  (32) 

RHS  of  (26)  =  2^1  +  ^£2  +  0(£3)^  (33) 

An  equation  for  b{T>  can  be  derived  by  substituting  results 
(32)  and  (33)  into  Eq.  (23)  for  k  —  1: 

e  +  5/32e3  +  0(ba)  e 

1  +  1  /4e2  +  0(b3)  ~  \+b<2>E2  +  0(e3) 

which  gives  us  b(2)  =  3 /32  and  thus,  we  obtain 

b=  ^7r(1  +  ^£2  +  °(£4))  '34') 
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It  is  desirable  to  express  the  normalized  Onsager  potential 
U (6)  as  a  function  of  the  normalized  polymer  concentra¬ 
tion  b.  For  that  purpose,  we  solve  for  e  from  (34)  and 
arrive  at 


s  —  ±, 


‘  b  —  (3/2)77 
(9/64)77 


U (6)  —  —  £  cos  20  —  —  cos  46  +  0{e 3) 


(35) 

(36) 


where  the  sign  ±  means  that  there  are  two  equivalent 
branches  of  nematic  state  near  b  =  (3/2)77.  These  two 
branches  are  equivalent  to  each  other  via  a  90  degree  rota¬ 
tion.  The  equilibrium  probability  density  is  given  by 


P(0)  = 


13  15 

- be  —  cos  26  +  e 2 - cos  46  +  0(e3) 

277  4b  64  b 


(37) 


where  e  is  given  in  (35).  This  is  the  first  nematic  branch  of 
polymer  ensembles  with  the  Onsager  potential.  Before  we 
plot  the  first  nematic  branch,  we  need  to  decide  on  how  to 
represent  the  nematic  state  as  a  function  of  b.  We  could  plot 
coefficients  (c2,  c4,  ...  ,  c2j, . . .)  as  functions  of  b.  But  that 
plot  is  a  curve  in  a  space  of  infinite  dimension,  which  is 
difficult  to  illustrate.  To  represent  the  phase  diagram  using 
a  simple  two-dimensional  plot,  we  plot  the  leading  non¬ 
zero  coefficient  in  (c2,  c4, ...,  c2j,  ■  ■  ■)■  This  approach  of 
plotting  the  leading  non-zero  Fourier  coefficients  is  in,  fact, 
consistent  with  what  has  been  done  in  the  case  of  nematic 
polymers  with  the  Maier-Saupe  interaction  potential.  For 
the  Maier-Saupe  potential,  c2  is  the  only  possible  non¬ 
zero  Fourier  coefficient  and  other  Fourier  coefficients  are 
always  zero  because  by  definition  the  Maier-Saupe  poten¬ 
tial  contains  only  the  cos  26  term.  In  contrast,  in  general  the 
Onsager  potential  contains  all  cos  2 jd  terms  where  j  goes 
from  1  to  infinity.  For  the  onsager  potential,  any  member  in 
set  (c2,  c4, ... ,  c2j, . . .)  can  be  the  leading  non-zero  coef¬ 
ficient.  For  the  first  nematic  branch  of  polymer  ensembles 
with  the  Onsager  potential,  c2  is  the  leading  non-zero  coef¬ 
ficient.  The  behavior  of  this  nematic  branch  is  shown  in 
Figure  1. 

In  Figure  1,  the  vertical  coordinate  c2  is  a  measure  of 
the  magnitude  of  the  Onsager  potential  for  the  first  nematic 
branch.  The  equilibrium  probability  density  is  related  to 
the  Onsager  potential  by  the  Boltzmann  distribution.  Thus, 
c2  also  measures  how  far  the  nematic  state  deviates  from 
the  isotropic  state.  In  this  sense,  c2  can  be  viewed  as  an 
unnormalized  order  parameter  for  the  first  nematic  branch. 
In  Figure  1,  the  solid  line  represents  a  very  accurate  numer¬ 
ical  solution,  which  is  treated  as  the  exact  solution,  and 
the  dashed  line  represents  the  asymptotic  solution  obtained 
above.  The  numerical  method  will  be  discussed  in  the  next 
section. 

For  polymer  ensembles  with  the  Maier-Saupe  interac¬ 
tion  potential,  the  polymer  orientation  distribution  is  often 
illustrated  by  plotting  in  polar  coordinate  system  a  func¬ 
tion  f{6)  showing  the  alignment  of  polymer  rods  along  the 


b 

Fig.  1.  The  first  nematic  branch  of  polymer  ensemble  with  the  Onsager 
interaction  potential.  This  nematic  branch  intersects  with  the  isotropic 
branch,  at  (3/2) it.  The  horizontal  coordinate  b  is  the  normalized  poly¬ 
mer  concentration;  the  vertical  coordinate  c2  is  the  coefficient  of  the 
cos  26  term  (the  leading  non-zero  term  for  the  first  nematic  branch)  in 
the  Fourier  expansion  of  the  normalized  Onsager  potential. 


direction  with  angle  6.  Function  f  (6)  is  defined  as  the  root- 
mean-square  of  the  inner-product  between  polymer  rods  in 
the  ensemble  and  the  direction  6. 


m  = 


J  J  ((cos  6,  sin  6)  ■  (cos  6,  sin  6))2p(6)dd 

cos  2{6  —  6)p{6)d6 


1 1 


1  +  cos  2(6  —  6) 


p(6)dd 


1  1  1 

2  +  2 


o2tt  „  „  „ 

cos  2 6  /  cos  2 6p(6)dd 
J  o 


11  1  377 

-  +  -a277cos20  =  ,/ -  4 - c2cos20  (38) 

2  2  2  V  2  8  b  2  V  ’ 


In  Figure  2,  we  plot  in  polar  coordinate  system  probabil¬ 
ity  density  p(6)  and  function  f(0)  of  a  typical  equilibrium 
state  on  the  first  nematic  branch.  Note  that  in  some  litera¬ 
ture,  the  graph  of  function  f{6)  was  mistakenly  identified 
as  an  ellipse.  This  misconception  was  probably  caused  by 
the  fact  that  the  function  form  of  f2(6).  when  written  in 
terms  of  the  Cartesian  coordinates  (x,  y),  has  an  elliptic 
form.  Using  x  —  cos  6  and  y  =  sin  6,  we  have 


f2(+y) 


1  377 

-  +  -c2cos2  6 


1  377  >  , 

-  +  — —  c2  )  cos‘  6  + 


8  b 
377 


377 


-  -  —  Cr 


1 

2  8  b 

377 

8 b' 


sin2 


6 

(39) 
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Fig.  2.  Probability  density  p(6)  and  function  f(6)  shown  in  polar  coor¬ 
dinate  system  for  a  typical  equilibrium  state  on  the  first  nematic  branch. 


It  can  be  shown  that  the  coefficient  (1/2—  (37t/(8/;))c2  is 
always  positive.  So  indeed  f2(x,y )  has  an  elliptic  form. 
But  the  graph  of  f(x,  y)  with  the  constraint  x2  +  y2  =  1  is 
not  exactly  an  ellipse  although  it  looks  not  very  different 
from  an  ellipse. 

In  Figure  2,  the  graph  of  probability  density  p(9)  is  sig¬ 
nificantly  different  from  an  ellipse.  It  has  a  dumbbell  shape. 
Actually,  the  graph  of  probability  density  p(8)  does  not 
always  have  a  dumbbell  shape.  For  small  value  of  c2,  that  is, 
for  b  sufficiently  close  to  (3/2)77,  the  graph  of  p(d)  looks 
like  an  ellipse  in  the  sense  that  it  is  convex  everywhere. 
In  particular,  it  is  convex  at  8  —  (1/2)7 7.  The  curvature  of 
p(0)  is  given  by 


«(0)  = 


p2(8)  +  2(p>(0)y--p(d)p"(0) 

(p2(d)  +  (p'(d)yf2 


(40) 


At  9  =  (1  /2)  77  and  for  small  value  of  c2  =  e,  we  have 


>’[l2,r)  =  t-te  +  0[el> 


1 


p'(-77)=0 


(41) 


2  77  j  p  ((1/2)77) 


(1  —  5e)  +  0(s2) 


Thus,  for  small  values  of  c2  =  e,  the  curvature  k  ((1/2)77) 
is  positive  and  the  graph  of  p(9)  is  convex.  The  results 
shown  in  Figure  1  demonstrate  that  the  asymptotic  solution 
is  still  valid  even  for  fairly  large  e.  Suppose  the  asymp¬ 
totic  expression  of  k  ((1/2)77)  given  in  (41)  is  valid  for 
moderately  large  s.  We  can  use  the  asymptotic  formula 
to  determine  at  what  value  of  e  the  curvature  changes 
sign  from  positive  to  negative.  Setting  k  ((1/2)77 )  =  0, 
we  obtain  e  =  0.2,  which  is  approximately  the  value  at 
which  p(8)  change  from  being  convex  to  being  concave 


at  6  —  (1/2)77.  This  has  been  confirmed  in  our  numerical 
simulations.  Once  p(6)  becomes  concave  at  8  =  (1/2)77, 
the  graph  of  p(9)  takes  the  shape  of  a  dumbbell. 


•  Case  2:  If  a2  =  0  but  a4  ^  0,  then  in  (29)  for  k  =  2  we 
must  have 


1- 


77 

2  M°) 


[42  -  1]  =  0 


which  yields  b (0)  =  (42— l)/277  =  (15/2)77.  Conse¬ 
quently,  for  k  >  2,  we  have  1  —  77/(2M°')[(2k)2  —  1]  = 
(42  —  (2k)1) /(A1  —  1)  <  0.  It  follows  from  (29)  that  a2k  —  0 
for  k  >  2.  As  we  did  in  case  1 ,  when  a4  ^  0,  we  re-define 
c4  =  a4e  +  j84e2  +  Off/)  as  the  new  e,  which  is  equiva¬ 
lent  to  setting  a4  =  1  and  /34  =  0.  Substituting  «4  =  I  and 
a2k  =  0  for  k  y  2  into  (25),  we  have 


RHS  of  (25) 


=  (s<x2k  +  e2p2k)TT  + 


k  =  4 

+  0(e2)  (42) 

otherwise 


As  in  case  1,  once  we  finish  calculating  coefficients  a2k, 
we  retain  both  the  first  order  terms  and  the  second  order 
terms  in  Eq.  (23)  to  calculate  coefficients  (32k .  We  first  cal¬ 
culate  b(l\  For  k  —  2,  Eq.  (23)  simply  can  be  written  as 

£+0(£3)  £+0(£3) 

l  +  0(£2)  ~  1  +  fiM1) 

which  immediately  implies  £>(1)  =  0.  For  k  =  4,  (23)  gives 
an  equation  for  /3g: 

£2/38  +  £2(l/4)  +  0(£3)  _  82-  1  £2/38 

l  +  0(fi2)  _  42  -  r  1  +  0(fi2) 

which  yields  /38  =  (5  /64).  For  k  {2,4},  Eq.  (23)  takes 
the  form 

£2fe  +  Q(£3)_  (2k)2 -1  E2P2k 
1  +  0(£2)  42  -  1  1  +  0(£2) 

which,  along  with  /34  =  0,  gives  us  (32k  =  0  for  k  y  4. 
Therefore,  the  nematic  state,  which  has  zero  coefficient  for 
cos  2 8  but  has  a  small  non-zero  coefficient  for  cos  48  in  the 
Fourier  expansion  of  U(9),  is  given  by 


U(8)  —  —  £  cos  40 - £2  cos  80+  0(e3) 

64 

b=  ^^(1  +  C>(£2)) 


(43) 


To  calculate  the  coefficient  of  the  second  order  term  in  the 
expansion  of  b,  we  retain  one  more  term  in  the  expansions 
of  (25)  and  (26) 


2 1  7 j 

RHS  of  (25)  =  £77  +  £3  — —  +  0(eA),  k  =  2 
128 

RHS  of  (26)  =  277  M  +  Ifi2  +  <2(£3) 


(44) 

(45) 
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An  equation  for  b(2>  can  be  found  by  substituting  results 
(44)  and  (45)  into  Eq.  (23)  for  k  =  2: 

s  +  (2I/128)e3  +  0(e4)  _  s 

1  +  (l/4)e2  +  0(fi3)  ~  1  +  b^e2  +  0{e3) 

which  gives  bl2)  —  11/128  and  therefore,  we  obtain 

l,=j’r(I+i^e2+0(‘4>)  <46> 


In  order  to  express  the  normalized  Onsager  potential  U (6) 
as  a  function  of  the  normalized  polymer  concentration  b, 
we  solve  for  s  from  (46)  and  arrive  at 


e  —  ±, 


'  &  —  (15/2)77- 
(165/256)7 r 


U(9)  =  — ecos4  9 - £2  cos80+0(£3) 

64 


(47) 

(48) 


Again,  the  sign  ±  indicates  that  there  are  two  equivalent 
branches  of  nematic  state  near  b  —  (15/2)7 t.  These  two 
branches  are  equivalent  to  each  other  via  a  90  degree  rota¬ 
tion.  The  equilibrium  probability  density  can  be  written  as 


p(6)  =  — +  £— cos40  +  £2  ^  £2cos80  +  (3(£3)  (49) 
277  4  b  256 b 

where  £  is  given  in  (47).  This  is  the  second  nematic  branch 
for  polymer  ensembles  with  the  Onsager  potential.  The 
behavior  of  this  nematic  branch  is  shown  in  Figure  3. 

In  Figure  3,  the  vertical  coordinate  c4  measures  the  mag¬ 
nitude  of  the  Onsager  potential  for  the  second  branch.  The 
equilibrium  probability  density  is  related  to  the  Onsager 
potential  by  the  Boltzmann  distribution.  Thus,  c4  also 


Fig.  3.  The  second  nematic  branch  of  polymer  ensembles  with  the 
Onsager  interaction  potential.  This  nematic  branch  intersects  with  the 
isotropic  branch  at  ( 1 5/2) 77.  Here  the  horizontal  coordinate  b  denotes 
the  normalized  polymer  concentration  whereas  the  vertical  coordinate 
c4  is  the  coefficient  of  the  cos  40  term  (the  leading  non-zero  term  for 
the  second  nematic  branch)  in  the  Fourier  expansion  of  the  normalized 
Onsager  potential. 


reflects  how  far  the  nematic  state  deviates  from  the 
isotropic  state.  In  this  sense,  c4  can  be  thought  of  as  an 
unnormalized  order  parameter  for  the  second  branch.  In 
Figure  3,  the  solid  line  corresponds  to  a  very  accurate 
numerical  solution,  which  is  taken  as  the  exact  solution,  and 
the  dashed  line  denotes  the  asymptotic  solution  obtained 
above.  Note  that  for  the  second  nematic  branch,  the  mean¬ 
ing  of  the  order  parameter  is  different  from  that  of  the  first 
nematic  branch.  For  the  first  nematic  branch  and  for  the 
case  of  the  Maier-Saupe  potential,  the  order  parameter  mea¬ 
sures  the  alignment  of  polymer  rods  along  one  direction 
(the  direction  along  which  the  alignment  is  maximized). 
For  the  second  nematic  polymer,  as  we  will  see  below,  the 
order  parameter  c4  measures  the  alignment  of  polymer  rods 
with  two  perpendicular  directions. 

Figure  4  depicts  in  polar  coordinate  system  the  probabil¬ 
ity  density  p{0)  of  a  typical  equilibrium  state  on  the  second 
nematic  branch.  Function  f(6),  the  root-mean-square  of  the 
inner-product  between  polymer  rods  and  the  direction  6,  is 
not  shown.  In  fact,  for  the  second  nematic  branch,  c2  —  0 
and  f(6)  =  const.  More  specifically, 


m  = 


1 1  377 

J  -  -| - c,  cos20  = 

V  2  8  b  2 


(50) 


Figure  4  shows  that  the  probability  density  p{9)  has 
four  fold  symmetry  and  the  polymer  rods  are  equally  likely 
being  aligned  along  the  horizontal  direction  and  being 
aligned  along  the  vertical  direction. 

•  Case  3:  If  a,  =  a4  =  0  but  a6  ^  0,  then  in  (29)  for  k  =  3 
we  must  have 


which  yields  M0)  =  ((62  —  l)/2)77  =  (35/2)77.  Conse¬ 
quently,  for  k  >  3,  we  have  1  —  (77/(2Z/0)))[(2A:)2  —  1]  = 
(62  —  (2k)2)/(62  —  1)  <  0.  It  follows  from  (29)  that  a2k  =  0 
for  k  >  3.  Since  a6  ^  0,  we  re-define  c6  =  a6s  +  /36s2  + 
0(s3)  as  the  new  s,  which  is  equivalent  to  setting  a6  =  1 


Fig.  4.  Probability  density  p(6)  shown  in  polar  coordinate  system  for 
a  typical  equilibrium  state  on  the  second  nematic  branch. 
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and  /36  =  0.  Substituting  a()  =  I  and  a2k  =  0  for  k  ^  3  into 
(25),  we  have 


RHS  of  (25) 

=  0«2  k  +  e2P2k)™  + 


k  =  6 

+  0(s3)  (51) 

otherwise 


Now  we  continue  to  calculate  coefficients  fi2k .  For  that  pur¬ 
pose,  we  keep  both  the  first  order  terms  and  the  second 
order  terms  in  Eq.  (23).  We  first  calculate  h(l>.  For  k  —  3, 
Eq.  (23)  simplifies  to 

e+0(e3)  £  +  0(s3) 

1  +  0(e2)  ~  1  +  eM1) 

which  gives  bil)  —  0.  For  k  —  6,  (23)  yields  an  equation  for 
Pl2- 


£2/312  +  £2(l/4)  +  <9(£3)  _  122-  1  £2/812 

l  +  0(£2)  ~~  62  -  1  ’  1  +  <9(e2) 

which  implies  /312  =  35/432.  For  k  {3,6},  Eq.  (23) 
becomes 


£2P2k  +  Q(E3)  _  (2k)2  -  1  e2/32 

1  +  0(£2) 


62  -  1  1  +  0(e2) 

which,  along  with  =  0,  yields  fi2k  =  0  for  k  ^  6.  There¬ 
fore,  the  nematic  state  that  has  zero  coefficients  for  cos  26 
and  for  cos  46  but  has  a  small  non-zero  coefficient  for 
cos  66  in  the  Fourier  expansion  of  U(6)  can  be  expressed  as 


35 

U  (6)  —  —  b  cos  66  —  yy2  cos  12  6+  0(s3) 
b=  y7T(l  +  0(£2)) 


(52) 


To  compute  the  coefficient  of  the  second  order  term  in  the 
expansion  of  b,  we  maintain  one  more  term  in  the  expan¬ 
sions  of  (25)  and  (26) 


14377 

RHS  of  (25)=  £  +  £3 — : - b<9(fi4), 

v  ’  864  v  ’ 


k  =  3 


RHS  of  (26)  =  277  (  1  +  -£2  +  <9(£3) 


(53) 

(54) 


An  equation  for  b(2>  can  be  derived  by  substituting  results 
(53)  and  (54)  into  Eq.  (23)  for  k  —  3: 

£+(143/864)£3  +  0(£4)  £ 


l  +  (l/4)£2  +  0(£3)  l  +  b<2)£2  +  0(£3) 

which  gives  us  0(2)  =  73/864.  So  putting  all  these  together, 
we  find 


b=  +  —  £2  +  0(£4) 

2  864  V  ’ 


(55) 


To  express  the  normalized  Onsager  potential  U(6)  as  a 
function  of  the  normalized  polymer  concentration  b ,  we 
solve  for  e  from  (55)  and  obtain 


£  =  ±, 


1  b-  (35/2)77 
(2555/1728)77 

35 


U(6)  —  — £cos60 - £  cos  120+  0(e  ) 

w  432 


(56) 

(57) 


The  sign  ±  indicates  that  there  are  two  equivalent  branches 
of  nematic  state  near  b  =  (35/2)77.  These  two  branches 
are  equivalent  to  each  other  via  a  90  degree  rotation.  The 
equilibrium  probability  density  is  given  by 


p(6)=  ^-  +  £^cos60  +  £2^^cos120+O(£3) 

(58) 

where  e  is  given  in  (56).  This  is  the  third  nematic  branch  for 
polymer  ensembles  with  the  Onsager  potential.  The  behav¬ 
ior  of  this  nematic  branch  is  shown  in  Figure  5. 

In  Figure  5,  the  vertical  coordinate  c6  gives  a  measure  of 
the  magnitude  of  the  Onsager  potential  for  the  third  branch 
while  the  equilibrium  probability  density  is  related  to  the 
Onsager  potential  by  the  Boltzmann  distribution.  Thus,  c6 
also  determines  how  far  the  nematic  state  deviates  from 
the  isotropic  state.  In  this  sense,  c6  can  be  treated  as  an 
unnormalized  order  parameter  for  the  third  branch.  As  we 
will  see  below,  for  the  third  branch,  the  order  parameter  c6 
actually  measures  the  alignment  of  polymer  rods  with  three 
directions.  In  Figure  5,  the  solid  line  denotes  a  very  accu¬ 
rate  numerical  solution,  which  is  treated  as  the  exact  solu¬ 
tion,  and  the  dashed  line  represents  the  asymptotic  solution 
obtained  above. 

In  Figure  6,  we  plot  in  polar  coordinate  system  probabil¬ 
ity  density  p(6)  of  a  typical  equilibrium  state  on  the  third 


Fig.  5.  The  third  nematic  branch  of  polymer  ensembles  with  the  Onsager 
interaction  potential.  This  nematic  branch  intersects  with  the  isotropic 
branch  at  (35/2)  it.  The  horizontal  coordinate  b  is  the  normalized  polymer 
concentration;  the  vertical  coordinate  c6  is  the  coefficient  of  the  cos  66 
term  (the  leading  non-zero  term  for  the  third  nematic  branch)  in  the 
Fourier  expansion  of  the  normalized  Onsager  potential. 
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Fig.  6.  Probability  density  p(6)  shown  in  polar  coordinate  system  for 
a  typical  equilibrium  state  on  the  third  nematic  branch. 


nematic  branch.  Function  f(6 ),  the  root-mean-square  of 
the  inner-product  between  polymer  rods  and  the  direction 
6,  is  not  shown.  For  nematic  branches  other  than  the  first 
nematic  branch,  we  have  c2  =  0  and  f(0)  —  const.  Figure  6 
demonstrates  that  the  probability  density  p(6)  has  six  fold 
symmetry  and  the  polymer  rods  are  equally  likely  being 
aligned  along  three  directions  that  are  uniformly  distributed 
on  the  unit  circle. 

The  process  of  finding  nematic  solutions  can  be  contin¬ 
ued  indefinitely.  Polymer  rod  ensemble  with  the  Onsager 
intermolecular  potential  has  infinitely  many  nematic 
branches.  The  first  nematic  branch  intersects  with  the 
isotropic  branch  at  b  =  (3/2)  77,  the  second  nematic  branch 
at  b  —  (15/2)7 t,  the  third  nematic  branch  at  b  —  (35/2)77. 
In  general,  if  a2  =  a4  —  ■  ■  ■  —  a2t/_i)  =  0  but  ct2j  ^  0,  then 
in  (29)  for  k  =  j  we  must  have 

which  yields  M0)  =  (((2/)2  —  l)/2)77.  Consequently, 
for  k  >  /,  we  have  1  —  tt  /  (2b<-0))[(2k)2  —  1]  = 
((2/)2  -  (2k)2)/((2jf  -  1)  <  0.  It  follows  from  (29)  that 
a2k  =  0  for  k  >  j.  Since  a2j  ^  0,  we  re-defme  c2j  = 
a2jE  +  (32jB2  +  0(e2)  as  the  new  £,  which  is  equivalent 
to  setting  a2j  =  1  and  /32,  =  0.  Substituting  a2,  =  1  and 
a2k  —  0  for  k  ^  /  into  (25),  we  have 


RHS  of  (25) 


=  ( sa2k  +  E1l32k)'JT  + 


k  =  2  j 


otherwise 


0(£3)  (59) 


which  gives  b (1*  =  0.  For  k  —  2j,  (23)  leads  to  an  equation 
for  p4j. 

£2P4j  +  e2(l/4)  +  0(e2)  _  (4  jf  -  1  e2p4j 

1  +  0(£2)  _  (2/)2  -  l'  1  +  0(£2) 

which  yields  p4]  =  (1/12)  ■  ((2y)2  -  l)/((2;)2).  For  k  £ 
{j,  2 j],  Eq.  (23)  has  the  form 

s2p2k  +  Q(s2)  _  (2k)2 -1  s2p2k 

l  +  0(£2)  (2  y')2  —  1  1  +  0(b2) 

which,  along  with  P2j  =  0,  implies  that  P2k  —  0  for  k  ^  2 j. 
Therefore,  the  nematic  state  that  has  zero  coefficients  for 
cos 26  through  cos 2 (j  —  1)0  but  has  a  small  non-zero  coef¬ 
ficient  for  cos  2  j6  in  the  Fourier  expansion  of  U  (6)  can  be 
described  mathematically  by 

1  (2  i)2  -l  . 

U  (6)  =  —£  cos  2/0—  —  ■  £-cos4,/0+O(£) 

b=  —————— 77  (l  +  0(£2))  (60) 

We  would  like  to  express  the  normalized  Onsager  potential 
U (6)  as  a  function  of  the  normalized  polymer  concentra¬ 
tion  b.  To  achieve  that  goal,  we  need  to  express  the  small 
parameter  £  in  terms  of  b.  Since  the  first  order  term  in  the 
expansion  of  b  with  respect  to  £  has  zero  coefficient,  to 
reveal  the  dependence  of  b  on  £,  we  need  to  find  the  leading 
term  with  non-zero  coefficient.  We  re-visit  the  expansions 
of  (25)  and  (26)  and  retain  one  more  term  in  the  expansions 

RHS  of  (25)  =  £  +  £3  (4^~  1}  77  +  0(£4),  k  =  j 

(61) 

RHS  of  (26)  =  277  M  +  -U2  +  <9(£3)  j  (62) 


Substituting  results  (61)  and  (62)  into  Eq.  (23)  for  k  =  j 
gives  us 

g  +  ((4(27')2—  l))/(24(2y)2)g3  +  0(g4)  = _ e _ 

l  +  (l/4)fi2  +  0(fi3)  l+M2»£2+0(£3) 

Solving  for  b{2)  from  the  equation  above  yields 


bd)  _  2(2y)~  + 1 
. .  24(2  j)2 


,  (2;)2 -  i  i  | 

b  —  - 77  1 


2(2 j)2  +  1  2  ,  ^  4, 


24  (2j)2 


£2  +  0(S) 


(63) 


Now  we  use  expansion  (63)  to  express  e  in  terms  of  the 
normalized  polymer  concentration  b  and  arrive  at 


To  calculate  coefficients  bil)  and  /j4/,  we  keep  both  the  first 
order  terms  and  the  second  order  terms  in  Eq.  (23).  For 
k  =  j,  Eq.  (23)  becomes 

£  +  0(fi3)  £+0(£3) 

l  +  0(£2)  ~  l  +  fi^d) 


e  —  ±. 


fe-(((2/)2-l)/2)77 


(((2/)2-1)(2(2/)2+1))/(48(2j)2)77 


(64) 


1  (2  i)2  -  1  , 

U(6)  =  -£  cos  2/0-  —  ■  (2/)2  £-cos4/0+O(£) 


(65) 
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Again,  the  sign  ±  implies  that  there  are  two  equivalent 
branches  of  nematic  state  near  b  —  (((2  j)2  —  l)/2)i7.  These 
two  branches  are  equivalent  to  each  other  via  a  90  degree 
rotation.  Once  the  Onsager  potential  is  calculated,  the  equi¬ 
librium  probability  density  is  given  by 


p(d)  =  —  +  £ — — - COS  2/0 

2-77  4b 


+  £' 


i(2j)2—  l)(4(2j)2—  1) 
48(27)2* 


cos4/0+O(e3)  (66) 


where  e  is  a  function  of  the  normalized  polymer  concen¬ 
tration  b  given  in  (64).  This  is  the  j-th  nematic  branch  of 
polymer  ensembles  with  the  Onsager  interaction  potential. 


4.  NUMERICAL  METHOD 

The  asymptotic  analysis  in  the  previous  section  tells  us 
two  things:  (i)  there  are  infinitely  many  nematic  branches 
for  polymer  ensembles  with  the  Onsager  interaction  poten¬ 
tial;  and  (ii)  the  j-th  nematic  branch  intersects  with  the 
isotropic  branch  at  b  —  (((2 j)2  —  1)/2)t7.  The  approximate 
expression  of  each  nematic  branch  near  the  intersection 
with  the  isotropic  branch  is  obtained  in  the  asymptotic  anal¬ 
ysis.  In  order  to  give  a  more  accurate  description  of  each 
nematic  branch  both  near  and  away  from  the  intersection 
with  the  isotropic  branch,  we  need  to  solve  (23)  for  an  arbi¬ 
trary  value  of  b.  In  Eq.  (23),  b  is  the  parameter  and  the 
infinite  sequence  of  coefficients  (c2,  c4, . . . ,  c2j, . . .)  is  the 
unknown.  From  physical  point  of  view,  it  is  desirable  to 
treat  (c2,  c4, ,  c2j, . . .)  as  functions  of  b.  In  the  phase 
diagrams  presented  in  the  previous  section,  for  practical 
reason,  we  show  only  the  leading  non-zero  Fourier  coeffi¬ 
cient  in  (c2,  c4, . ..,  c2j, . . .)  as  a  function  of  b.  Specifically, 
suppose 

c2(b)  =  •  •  •  =  c2(M)(fc)  =  0,  c2j(b)  =  r(b) 

In  phase  diagrams,  we  plot  r(b )  as  a  function  of  b  and 
indicate  that  r(b)  is  the  coefficient  for  cos 2/0.  Mathemat¬ 
ically,  however,  it  is  difficult  to  solve  for  r(b)  in  Eq.  (23) 
for  a  given  value  of  b.  There  are  two  mathematical  dif¬ 
ficulties:  (i)  as  in  the  case  of  nematic  polymers  with  the 
Maier-Saupe  potential,  r(b)  is  not  a  single-valued  func¬ 
tion  of  b  and;  (ii)  function  r(b)  is  not  defined  at  all  for 
b  <  (377/2),  that  is,  nematic  state  does  not  exist  at  all 
before  the  first  nematic  branch  appears.  Practically,  for  each 
nematic  branch,  it  is  more  convenient  to  treat  r  as  the  inde¬ 
pendent  variable  and  treat  b  as  a  function  of  r  instead. 
This  convenient  mathematical  formulation  has  also  been 
adopted  for  the  study  of  nematic  polymers  with  the  Maier- 
Saupe  potential.3-45^7  For  the  case  of  the  Maier-Saupe 
potential,  the  value  of  r  completely  determines  an  equi¬ 
librium  state.  We  believe  that  for  the  case  of  the  Onsager 
potential,  this  is  still  true:  the  value  of  r  completely  deter¬ 
mines  an  equilibrium  state  on  each  given  nematic  branch. 


This  assertion  has  been  confirmed  in  our  numerical  simu¬ 
lations.  But  a  rigorous  mathematical  proof  for  the  assertion 
is  still  open. 

Now  we  present  the  numerical  method  we  use  to  solve 
non-linear  system  (23).  For  the  convenience  of  presenta¬ 
tion,  we  introduce  the  vector  notation 


C  =  (c2,  c4, ... ,  c2k, . . .)  (67) 


Below  we  describe  an  iterative  method  for  calculating  the 
first  nematic  branch.  Similar  procedure  can  be  employed  to 
calculate  other  nematic  branches.  As  we  discussed  above, 
in  our  new  mathematical  formulation,  we  treat  c2  —  r  as 
the  independent  variable  and  treat  b  and  (c4, . . . ,  c2k, . . .) 
as  functions  of  r.  In  the  numerical  method,  r  is  the  input, 
and  b  and  (c4, . . . ,  c2k, . . .)  are  the  output.  More  precisely, 
we  start  with  a  given  value  of  r  and  proceed  to  calculate  the 
corresponding  values  of  b(r)  and  (c4(r), . . . ,  c2k(r), . . .). 
Here  is  an  outline  of  the  iterative  method: 

•  Step  0:  We  set  n  =  0  and  set  C(0)  =  (r,  0, . . . ,  0, . . .). 
Here  n  is  the  counter  recording  the  number  of  iterations 
completed  so  far  and  C,ri)  is  the  starting  vector  for  iteration. 

•  Step  1:  We  use  C<n>  —  (c2"\  c4"\  . . . ,  c2k\  . . .)  to 
calculate 


,„)  _  f02n  cos2A-0exp[E”  i  4°  cos2 id]dd 
fT  expE^i  cos2 i0]dO 

If  Q2k  is  the  solution  of  Eq.  (23),  then  for  k  =  1,  Eq.  (23) 
becomes 

Q<n)  =  —  r 

Recall  that  in  our  mathematical  formulation,  r  (i.e.,  coef¬ 
ficient  c2)  is  the  independent  variable  and  b  is  a  function 
of  We  set  the  value  of  b  to  satisfy  this  equation. 

•  Step  2:  We  calculate  the  new  value  of  b 


b(n)  = 


Sir 

4  Q?r 


•  Step  3:  For  k  >  1,  we  use  Eq.  (23)  with  the  values  of 
Q2k  and  to  update  coefficient  c2k: 


(n+\) 

C2A 


4^2A-  b(„) 

(4k2  —  1)77 


for  k  >  1 


•  Step  4:  We  update  the  coefficient  vector. 


C("+1)  =  (r>  c4"+1), . . . ,  c2£+1),  . . .) 


•  Step  5:  We  use  |  |C(”+1)  —  C(,,)  1 1  to  determine  if  the  itera¬ 
tion  has  converged.  If  |  |C/"+1)  —  C(,,)  1 1  is  less  than  the  given 
tolerance,  then  we  stop  the  iteration  and  output  the  most 
recent  coefficient  vector  C(”+1)  and  b ^  as  the  solution  of 
Eq.  (23).  Otherwise,  we  increase  the  counter  by  one  (set 
n  —  n+  1)  and  repeat  the  iteration  starting  at  Step  1. 
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Fig.  7.  Comparison  of  the  first  nematic  branch  for  the  Onsager  poten¬ 
tial  and  the  only  nematic  branch  for  the  Maier-Saupe  potential. 


Because  the  Smoluchowski  equation  is  a  diffusion  equa¬ 
tion  we  expect  the  probability  density  and  the  Onsager 
potential  to  be  smooth  functions.  For  a  smooth  function, 
the  Fourier  series  expansion  converges  exponentially  as 
the  number  of  terms  included  in  numerical  solutions  is 
increased.  Therefore,  numerical  solutions  accurate  up  to 
the  computer  precision  can  be  obtained  using  the  numer¬ 
ical  method  described  above  with  only  moderately  large 
number  of  Fourier  terms  included.  The  numerical  results 
are  shown  in  Figures  1-6.  Figure  7  compares  the  phase 
diagram  for  the  Onsager  potential  (showing  only  the  first 
nematic  branch)  and  the  phase  diagram  for  the  Maier- 
Saupe  potential.  As  shown  in  Figure  7,  the  Maier-Saupe 
potential  predicts  that  the  isotropic  to  nematic  phase  tran¬ 
sition  occurs  at  b  =  4  whereas  the  Onsager  potential  pre¬ 
dicts  that  the  isotropic  to  nematic  phase  transition  occurs 
at  b  =  (3/2)77^4.7124. 


5.  FREE  ENERGY  CALCULATION  AND 
STABILITY  OF  THE  ISOTROPIC  BRANCH 
AND  NEMATIC  BRANCHES 

In  the  Doi-Hess  model,  the  polymer  orientation  distribu¬ 
tion  is  described  by  the  probability  density  p(6).  The  free 
energy  of  polymer  orientation  distribution  p(9)  is 

Free  energy  of  orientation  distribution  p(6) 


G([P ])  = 


kRT 


=  fjpW  (\ogp(d)+l-U(6))dd 


(68) 


where  U (6)  is  the  normalized  Onsage  interaction  potential 

U(6)  =  b  \sm(6  -  6)\p(0)  dd  (69) 

Jo 


For  the  isotropic  branch,  the  Onsager  potential  is  U (6)  —  0, 
the  probability  density  is  p(6)  —  1/(2 tt)  and  the  corre¬ 
sponding  free  energy  is 


^Isotropic  l0g(277) 


(70) 


Note  that  in  the  expression  for  the  Onsager  potential 
U(6)  =0  we  have  dropped  the  constant  term  4b/ (2ir). 
If  we  keep  this  constant  term,  then  the  free  energy  of 
the  isotropic  state  is  GIsotropic  =  ( b/TT )  —  log(2ir),  which 
increases  linearly  with  the  normalized  polymer  concentra¬ 
tion  b.  Since  the  stability  of  an  equilibrium  state  is  deter¬ 
mined  by  the  free  energy  landscape  at  a  fixed  value  of  b, 
dropping  a  /^-dependent  constant  term  will  not  change  the 
stability. 

For  an  anisotropic  orientation  distribution  (not  necessar¬ 
ily  an  equilibrium  state),  the  Fourier  series  expansion  of 
probability  density  p(6)  and  that  of  the  normalized  Onsager 
potential  U (6)  are  related  by 

P(0)=  ^-  +  Efl2yCos2/0 
217  7=1 

oo 

U(0)  =  —  ^2  c2j  cos  2jd  (71) 

7=1 

4b 

cv  =  avj2]y^i 

Here  it  is  important  to  point  out  that  relation  (71)  is  true 
for  all  probability  densities.  In  particular,  the  probability 
density  p(6)  does  not  have  to  be  an  equilibrium  state.  If 
p(6)  is  an  equilibrium  probability  density,  then  in  addition 
to  relation  (71),  we  also  have  the  Boltzmann  relation 


log  P(0)  =  -  U(6)  -  log(Z)  =  J2  c2 j cos  2 jd  -  log (Z) 


7=1 


where  the  partition  function  Z  is  given  by 

Z  =  J  exp  c2j  cos 2 jd'j  dO 

The  free  energy  of  an  equilibrium  state  is 

c271  1 1  00 

G([p])  =  Jo  p(d)l -J2c2j cos2j 6- log(Z) 

7T  °° 

=  7rEC27fl2y-log(Z) 


(72) 

(73) 


d6 


7=1 


=  ^E4j((2j)2-l)-lo  g(Z)  (74) 


7=1 


For  the  first  nematic  branch,  the  Onsager  potential  has  the 

form  U{6)  =  —  (ecos2P+  (l/16)e2cos40H - )  and  the 

partition  function  is 


dd 


r2ir  /  1 

Z=  /  exp  I  e  cos  20 -| - £2cos40H - 

J  o  V  16 

=  / 

Jo 


2ir  /  1  11 

1  +  -£2  cos2  26+  - - £4  cos2  46 

2  2  162 


1  3 

-  •  — £4  cos2  26  cos  46 
6  16 
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24 


s4  cos4  20  +  ■ 


dO 


/  1  ,  25 

=  2tc  I  1  H — e2  H - e4  +  •  ■ 

v  4  1024 

The  log  partition  function  takes  the  form 
log(Z)  =  log(277)  +  log  (  1  +  ^E2  + 


=  l0g(277)- 


—  £2 
4 


25 

1024* 


1  /I 


2  l  4S 


=  log(277 )  +  I  —  £2 - £4 

S  \4  1024 


(75) 


+  ••• 
(76) 


The  free  energy  of  an  equilibrium  state  on  the  first  nematic 
branch  is 


Fig.  8.  Free  energy  of  equilibrium  state  on  the  first  nematic  branch  as 
a  function  of  the  normalized  polymer  concentration  b. 


IT 


15 


;  =  —  3£2  +  — £4  -  —  £2 - 1 

Fust  8  b\  256  /  \4  1024 


1 


-l0g(277)- 


3tt 
8 b 


£2  H - - —  £4  — 

{z1-— e4] 

^  +  -1 

|_  256 

^  256  ) 

377 

-l0g(277) 


(77) 


Using  the  fact  that  for  the  first  nematic  branch  b  —  (377)/ 
2(1  +  (3/32)£2),  we  obtain 


r  =  — 

F”st  8  b 


e2  + 


5 

256 


s4-le2  + 


17 

256* 


-l0g(277) 


377  12  , 

=  8/?  256S  -l0g(2,,)+  ^ 

Q  77  ] 

= - £4  —  log(277 )  H 

512  b  ’ 


977  1  /  b  —  (3/2)77 

512  ~b  \  (9/64)77 

4"  •  •  •  <  ^Isotropic 


-l0g(277) 


(78) 


Since  GIsotropic  +  log(277)  =  0,  for  simplicity  of  presenta¬ 
tion,  we  plot  G  +  log(277 )  as  the  free  energy.  Figure  8 
shows  the  free  energy  of  the  first  nematic  branch  as  a 
function  of  the  normalized  polymer  concentration  b.  The 
solid  line  represents  a  very  accurate  numerical  solution 
obtained  using  the  numerical  method  described  in  the  pre¬ 
vious  section,  which  is  treated  as  the  exact  solution,  and 
the  dashed  line  represents  the  asymptotic  solution  obtained 
above.  It  is  clear  that  the  asymptotic  solution  is  a  very 
good  approximation  to  the  exact  solution  over  a  very  large 
range  of  b.  Actually,  the  difference  between  the  asymptotic 
solution  and  the  exact  solution  will  not  increase  when  b  is 
increased  further  (see  Fig.  9). 

In  Figure  8,  the  dotted  line  shows  the  free  energy  of  the 
isotropic  branch.  Notice  that  the  free  energy  of  the  first 
nematic  branch  is  below  the  free  energy  of  the  isotropic 


branch.  The  first  nematic  branch  appears  at  b  =  (3/2)77. 
For  b  >  (3/2)77  but  close  to  (3/2)77,  the  first  nematic 
branch  is  the  only  nematic  branch.  Therefore,  the  lower 
free  ehergy  of  the  first  nematic  branch  implies  stability.  For 
b  <  (3/2)77,  the  isotropic  branch  is  the  only  solution.  As  a 
result,  for  b  <  (3/2)77,  the  isotropic  branch  is  stable.  The 
stability  of  the  isotropic  branch  for  b  >  (3/2)77  will  be  ana¬ 
lyzed  later  in  this  section. 

For  the  /-th  nematic  branch,  the  Onsager  poten¬ 
tial  has  the  form  U(6 )  =  —  (ficos2/0  +  ((2/)2  —  1)/ 
( 1 2 (2/)2 )fi2 cos4/Yi  +  •  ■  ■)  and  the  partition  function  is 


(2/)2  -  1 


Z  =  /  exp  (  £  cos  2/0  + 

4>  \  12(2y)2 

-fH 


£2  cos 4/0+  ■ 


dd 


—  £2  cos2  2/0  +  -  •  {{2j)~  ^  £4 cos2  4 jd 

2  J  2  122(2y)4  J 


3  (2y)2 —  1 
6 

1 


12(2y)2 


£4  cos2  2/0  cos  4/0 


+  —  £4  cos4  2/0  H — 


dd 


=  277  |  1  +  -^£ 


,2  ,  (4(2/)2-l)2e4  + 


9-64(2  /)4 

The  log  partition  function  has  the  expression 

,  ^  ,  1  2  ,  (4(2/)2- 1)2 

log(Z)  =  log(277)+log^l+4£  +  9.64(2/)4 

,  ,  ,  (l  2  ,  (4(2/)2  —  l)2  4 

’  ^4  9-64(2 /)4 

1/1 


(79) 


£4  +  *  *  * 


2\4E 


=  l0g(277) 


lfi2  2((2/)2  +  2)2-9  4 

4  9-64(2 /)4  1 


(80) 
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Fig.  9.  Free  energies  of  the  first  three  nematic  branches  as  functions 
of  the  normalized  polymer  concentration  b.  The  dotted  line  is  the  free 
energy  of  the  isotropic  branch. 


The  free  energy  of  an  equilibrium  state  on  the  first  nematic 
branch  is  given  by 


7 7 

Gj-th  -  8 b 


((2;)2-l)£2  + 

I 


(4(2y)2—  l)((2y)2—  l)2 
9-16(2;)4 


.  2  2((2y)2  +  2)z  — 9  41  ,  , 

A  9-64(2  ;)4  '  ’ 


((2y)2  —  1)tt 
8  b 


c2  !  4(2y)4  —  5(2j)2  + 1  „4 


-  £  - 


9  -  16(2j')4 

2((2;)2  +  2)2  — 9 


9-16(2;)4 


2b 


((2y)2  —  1)tt- 


-l0g(277) 


(81) 


Using  the  fact  that  for  the  ;-th  nematic  branch 
b  —  ((2;)2  —  l)/27r(l  +  (2(2;)2  +  1)/(24(2;)V  + 

(9(fi4)),  we  obtain 


Gj.  *  = 


((2y)2  —  1)tt 
8  b 


4(2y)4  —  5(2;)2  +  1 
9-16(2;)4 

2  2((2y')2  +  2)2  —  9 


916(2;) 


A4 


2b 


((2y)2  —  1)77- 


-  log(2tr) 


((2;)2-l)77  6(2;)2  +  3  4  ,  , 

8  b  9-16(2  j)2  J 

2(2  j)4  —  (2;)2  —  1  1  , 

- - - - — - - - 77-  —  £4  —  log  (2  77 )  4~ 

3  ■  128(2;)2  b  eV  ’ 


2(2;')4  —  (2;)2  —  1 
384(2 j)2 


1/  b  —  (((2;)2  —  l)/2)77  y 

b  V (((2;)2  —  1)(2(2;)2  +  1))/(48(2;)2)77 / 

-  l0g(277)  +  -  -  - 

^Isotropic  (82) 

Figure  9  shows  the  free  energy  of  the  first  three  nematic 
branches.  It  is  clear  that  for  all  nematic  branches,  the  free 
energy  is  below  that  of  the  isotropic  branch,  which  was 
established  analytically  in  (82)  above.  The  first  nematic 
branch  has  the  lowest  free  energy.  The  higher  the  index 
the  higher  the  free  energy.  The  free  energy  of  the  second 
nematic  branch  is  above  that  of  the  first  branch  and  the 
free  energy  of  the  third  nematic  branch  is  above  that  of  the 
second  branch  and  so  on. 

Figure  9  also  demonstrates  that  the  asymptotic  expres¬ 
sion  (82)  is  valid  for  a  very  large  range  of  b  and  valid  for 
all  nematic  branches.  The  asymptotic  solution  is  derived 
based  on  the  assumption  that  the  amplitude  of  the  Onsager 
potential  is  small.  So  it  is  still  a  mystery  why  the  differ¬ 
ence  between  the  asymptotic  solution  and  the  exact  solution 
remains  small  as  the  amplitude  of  the  Onsager  potential 
is  increased  over  a  very  large  range  corresponding  to  the 
increase  in  b.  The  mathematical  explanation  behind  this 
mystery  will  be  investigated  in  a  future  study. 

We  use  free  energy  calculation  to  study  the  stability  of 
the  isotropic  branch  for  b  >  (3/2)77  (i.e.,  after  the  first 
nematic  branch  appears).  The  free  energy  diagram  (Fig.  9) 
suggests  that  the  isotropic  state  is  unstable  for  b  >  (3/2)77. 
To  show  the  instability,  it  is  sufficient  to  show  that  it  is 
unstable  with  respect  to  a  perturbation  of  one  particular 
form.  We  consider  a  perturbation  to  the  isotropic  state  of 
the  form 

P(0)=  ^-(l  +  i?  cos  20)  (83) 

277 

where  17  is  a  small  parameter.  The  perturbed  probability 
density  satisfies  /0°°  p(d)dO  =  1.  Since  in  general  the  per¬ 
turbed  probability  density  is  no  longer  an  equilibrium  state, 
we  cannot  apply  the  Boltzmann  relation.  To  calculate  the 
free  energy  of  the  perturbed  probability  density,  we  can¬ 
not  use  the  expression  given  in  (74)  for  equilibrium  states. 
Instead  we  have  to  go  back  to  the  free  energy  expression 
given  in  (68)  for  a  general  polymer  orientation  distribu¬ 
tion  (not  necessarily  an  equilibrium  state).  Recall  that  the 
Fourier  series  expansion  of  p{0)  and  the  Fourier  series 
expansion  of  U(6)  are  connected  by  (71).  The  Onsager 
potential  corresponding  to  the  perturbed  probability  den¬ 
sity  is 

(yi  \  4/?  2/7 

—  )• — cos  20  — —V — cos  2 6  (84) 

277/  3  377 

The  second  component  (interaction  potential)  of  the  free 
energy  expression  given  in  (68)  is 
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fo  p(0)(^U(6)\de  =  p{0)(-r,^cos2o\d9 


71  b  ,  b 

=  -—■71—  77=7)^—  (85) 

277  377  077 

The  first  component  (entropic  contribution)  of  the  free 
energy  expression  given  in  (68)  is 


n2.1T 

/  p(9)\ogp(6)d0 

Jo 

n2ir 

—  /  p(6)\og(l  +  7]cos26)d9  —  log(27r) 

J  o 

—  J  p(9)  ^T7Cos2 9—  ^-cos220  +  O(t]3)\  dd  +  G 


=  ■  71- 77-  \  \  -l0g(277)+  0(173) 

277  2  2 

772 

—  +  0(v 3)  +  Gisotropic 


isotropic 


(86) 


The  total  free  energy  of  the  perturbed  probability  density  is 


We  expand  the  free  energy  as  a  power  series  of  rj  and  keep 
only  the  constant  terms  and  172  terms.  Note  that  all  17  terms 
sum  to  zero  because  the  unperturbed  probability  density 
is  an  equilibrium.  The  expansion  of  the  first  component 
(entropic  contribution)  of  free  energy  is 


(p{e)(9)  +  7] cos 29)  log  (p(e>(9)  +  i7cos20)  d9 


=  fj  (pw(0)logPw(0) 

+  ^^(0)(T,COS20)2  +  O(7?3))^ 

n2.1T 

=  /  p(e)(0)logpw(0)</0 

j 0 

1  r277  1 

+  4VJ0  ^y(!  +  cos4^+0(p3)  (90) 


The  Onsager  potential  corresponding  to  the  perturbed  prob¬ 
ability  density  is 


G(t?)  =  ^isotropic  2  )  +  G(r,3)  (87) 

For  fixed  b  >  (3tt)/2  and  small  17,  we  have  G(p)  < 
Gisotropic ’  which  implies  that  for  b  >  (3ir)/2  the  isotropic 
state  is  unstable. 

So  far  we  have  concluded  that  (i)  the  isotropic  state  is  the 
only  solution  for  b  <  (3ir)/2  and  is  stable  for  b  <  (3t7)/2; 
(ii)  the  first  nematic  branch  appears  at  b  —  (3ir)/2  and  is 
stable  for  b  >  (3tt)/2.  The  y-th  nematic  branch  appears  at 
b  =  (((2y)2—  1) /2)77.  Although  the  free  energy  diagram 
(Fig.  9)  above  does  not  yield  a  definite  conclusion  on  the 
stability  of  the  y-th  nematic  branch  for  j  >  1,  intuitively 
it  is  reasonable  to  expect  that  for  j  >  1  the  y-th  nematic 
branch  is  unstable.  Below  we  use  free  energy  calculation 
to  show  the  instability.  Again,  to  conclude  the  instability, 
we  only  need  to  identify  one  unstable  model.  For  j  >  1, 
the  unperturbed  equilibrium  probability  density  on  the  y-th 
nematic  branch  has  the  form 

1  00 

p{e)(9)=—  +  Y.a2kcos2k9  (88) 


U(9)  —  U^e)  (0)  —  17  •  —  cos  20  (91) 

where  U(e\9)  is  the  Onsager  potential  corresponding  to  the 
unperturbed  equilibrium  probability  density  given  by 

00  4  h 

t/W(0)  =  -^>  COS  2^0  (92) 

k= 2  1 

We  expand  the  second  component  (the  Onsager  potential) 
of  the  free  energy. 

1  r2lT  (  4  b  \ 

-  J  (p(c)  (0)  +  7]  cos  20)  (  G(f)  (0)  -  v  ■  y  cos  20  j  d0 

=  -  J  l  p{e)(9)U{e)(0)d9-— 7]2  cos2  29  jd0 

=  \  r  P(e) ( 0) G(e) (9)d9—  -^-7]2  (93) 

2  Jo  3 

Adding  together  the  two  components  of  the  free  energy,  we 
have 


Specifically,  the  unperturbed  probability  density  satisfies 
a2  =  0.  We  consider  a  perturbation  to  the  equilibrium  prob¬ 
ability  density  of  the  form 

p(0)  =  p(e)(0)  +  T7cos20  (89) 

where  17  is  a  small  parameter.  The  perturbed  probability 
density  satisfies  f'77  p (0) d 9  —  1.  To  proceed  with  the  free 
energy  calculation,  we  first  look  at  the  expansion  of  x  ■  log  x 


G(i?)  =  G(0)-i72772 

/  2b  1  r2lr  1  \ 

7r-l  ^a+c»«v«)+0(,-') 

(94) 

For  the  y-th  nematic  branch,  the  equilibrium  probability 
density  is 


(*  +  <5)  Tog(v  +  5)  =  x  •  log  x  +  (log  x  +  1)5-| - 52H - 

2  x 


p(e)  =  —  (1  +  sq2j  cos2jd  +  e2r2j  cos 4/0  +  0(s3))  (95) 
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where  b  and  coefficients  q2j  and  r2;  are  given  by 
0-(((2;)2-1)/2)t7 


£  =  d=_ 


(((2/)2  —  l)(2(2j)2  +  1))/ (48(2;)2)t7 


q2j  =  2  77 


(2y)2 -  1 


(96) 


rV  = l7T 


1  "  4  b 

((2j)2  —  l)(4(2j)2  —  1) 
48(2  ;)20 


The  reciprocal  of  the  equilibrium  probability  density  has 
the  expansion. 

(e)^  =  2ir(l  —  eq2j  cos2;0  —  e2r2j  cos4j6 

+  e2qlj  cos2  2j6+  0(e3))  (97) 

We  calculate  the  second  term  in  the  coefficient  of  rj1  in 
(94). 


1  r2”  l 

— -  - - -(l+cos4  6)d0 

477“  J 0  p(c>(0) 

1  +  e2(^2-^2)+"-  for;=l 

=  '  l-e^4  +  «2^4  +  -"  for;  =  2  (98) 


[  1  +£2^2;  +  "'  for;  >2 

The  ;-th  nematic  exists  only  for  b>  (((2;)2-1)/2)t7. 
Thus,  for  ;  >  1  and  for  b  near  (((2;)2  —  l)/2)77  (i.e.,  for 
nematic  states  of  small  amplitude  e),  we  have 

G(v)  <  G( 0)  -  7]~ 7r~  — -  -  1  +  0(e)j  +  0( 173) 

<  G(0)  -  2t/2tt2  +  0(p3)  <  G(0)  (99) 

which  implies  that  for  ;  >  1,  the  ;-th  nematic  branch  is 
unstable.  In  the  above,  we  have  concluded  the  stability  of 
the  first  nematic  branch  based  on  the  fact  that  it  is  the  equi¬ 
librium  state  with  the  lowest  free  energy.  Nevertheless,  it 
is  worthwhile  to  examine  directly  the  stability  of  the  first 
nematic  branch  in  the  presence  of  a  small  perturbation.  We 
consider  perturbed  probability  densities  of  the  form 


p(0)  =  pw(0)  +  A2ltcos2A0  (100) 

k=  1 

It  can  be  shown  that  including  Fourier  modes  cos(2A-f- 1)0 
and/or  sin  k6  in  the  perturbation  will  only  increase  the 
free  energy.  Consequently,  to  derive  the  stability  we  only 
need  to  consider  perturbations  consisting  of  Fourier  modes 
cos  2 kO.  The  Onsager  potential  corresponding  to  the  per¬ 
turbed  probability  density  is 

00  4b 

£/(fl)  =  t/w(e)-7?£  ,  A2t. cos 2k 6  (101) 

k=  l  (2*)2-l 


where  U(e\0)  is  the  Onsager  potential  corresponding  to  the 
unperturbed  probability  density.  We  expand  the  two  com¬ 
ponents  of  the  free  energy  as  a  power  series  of  r/.  Again, 
we  collect  only  the  constant  terms  and  172  terms.  The  sum 
of  all  17  terms  must  be  zero  because  the  unperturbed  prob¬ 
ability  density  is  an  equilibrium  state. 
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We  calculate  7', ,  7),  7) ,  and  74  one  by  one.  7’,  is  bounded  by 
r_  26”  A2,  >  2b-\2k_  2b  ~ 

1  77  “  (2A)2  -  1  "  77  £  15  1577  £  24 

(103) 


For  b  near  (3tt)/2  (i.e.,  for  nematic  states  of  small  ampli¬ 
tude  b  on  the  first  branch),  we  have 
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00 
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2k 


For  b  near  (37r)/2,  we  use  (0  —  (3tt)/2)  =  s2(9tt)/64  and 
result  (98)  to  expand  7)  with  respect  to  b 

,  /  1  /■2,r  cos2  20  20 

7)  =  A,  ( — -  — - r/0 - 

"  y  277“  7o  P(C  (0)  377 


(104) 


At  0  =  (3tt)/2  where  the  first  nematic  branch  emerges 
from  the  isotropic  branch,  coefficients  q2  and  r2  take  the 
values  q2=  1  and  r2  =  5/16,  which  leads  to 
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It  follows  that  for  b  near  (37 t)/2,  T2  has  a  lower  bound 
given  by 

T2  >  ^£2A2 

For  b  near  (3t7)/2,  the  nematic  equilibrium  probability 
density  on  the  first  branch  has  the  expansion  given  in  (95). 
Using  that  expansion,  we  derive  lower  bounds  for  7)  and  T4. 
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Substituting  all  these  results  into  (102),  we  arrive  at 
G(V)>  G(0)  +  tj2^ 


4£a2,-V£2a2.  £a2,+3£2a: 


>  G(0)  +  i,2yl2^A^£2A2  >G(0) 


(106) 


which  concludes  the  stability  of  the  first  nematic  branch. 


6.  CONCLUSIONS 

We  have  carried  out  an  asymptotic  analysis  on  the  phase 
diagram  of  nematic  polymer  monolayers  with  the  Onsager 
interaction  potential.  In  the  case  of  polymers  with  the 
Maier-Saupe  potential  which  is  an  approximation  to  the 
Onsager  potential,  there  is  only  one  nematic  branch.  In 
contrast,  polymers  with  the  Onsager  interaction  poten¬ 
tial  have  infinitely  many  nematic  branches.  For  small 
b  (the  normalized  polymer  concentration),  the  isotropic 
state  is  the  only  equilibrium.  As  b  is  increased,  these 
nematic  branches  appear  one  by  one,  starting  with  the  first 
nematic  branch  at  b  —  (3/2)77,  then  the  second  nematic 
branch  at  b=  (15/2)77,  then  the  third  nematic  branch 
at  b  —  (35/2)77.  The  /-th  nematic  branch  appears  at 
b  =  ((2 j)2  —  For  nematic  states  on  the  first  branch, 


polymer  rods  are  aligned  along  one  direction,  similar  to 
nematic  states  in  the  case  of  the  Maier-Saupe  potential. 
Indeed,  the  phase  diagram  we  obtained  shows  that  the  first 
nematic  branch  is  close  to  the  only  nematic  branch  for  the 
Maier-Saupe  potential.  For  nematic  states  on  the  second 
branch,  the  behavior  of  the  polymer  orientation  distribu¬ 
tion  is  drastically  different.  Instead  of  being  aligned  along 
one  direction,  polymer  rods  are  aligned  with  two  perpen¬ 
dicular  directions  and  the  probability  density  has  four  fold 
symmetry.  For  nematic  states  on  the  /-th  branch,  polymer 
rods  are  aligned  with  /  directions  uniformly  distributed 
on  the  unit  circle,  and  the  probability  density  has  2 /  fold 
symmetry.  We  use  the  leading  non-zero  coefficient  in  the 
Fourier  expansion  of  the  Onsager  potential  as  an  unnor¬ 
malized  order  parameter.  For  the  first  nematic  branch,  the 
order  parameter  measures  the  alignment  along  the  director 
(the  direction  along  which  the  alignment  is  maximized), 
similar  to  the  case  of  the  Maier-Saupe  potential.  For  the 
subsequent  nematic  branches,  however,  the  order  parameter 
should  be  viewed  as  a  measure  on  the  deviation  of  the  prob¬ 
ability  density  from  the  isotropic  distribution.  We  have  fur¬ 
ther  designed  a  spectrum  method  for  calculating  the  phase 
diagram  of  nematic  polymers  with  the  Onsager  potential. 
The  spectrum  method  is  capable  of  producing  numeri¬ 
cal  solutions  that  are  accurate  up  to  the  computer  preci¬ 
sion,  which  we  use  as  the  exact  solution.  All  asymptotic 
results  obtained  in  this  paper  are  confirmed  by  the  numer¬ 
ical  simulations.  Although  mathematically  the  asymptotic 
analysis  applies  only  to  nematic  states  with  the  Onsager 
potential  of  small  amplitude,  the  comparison  of  the  asymp¬ 
totic  results  and  numerical  results  demonstrates  that  the 
asymptotic  expressions  derived  are  still  valid  for  a  very 
large  range  of  b,  even  when  the  amplitude  of  the  Onsager 
potential  is  very  large.  We  have  also  performed  free  energy 
calculations  to  study  the  stability  of  the  isotropic  branch 
and  the  nematic  branches.  Before  the  appearance  of  the 
first  nematic  branch  at  b  —  (3/2)77,  the  isotropic  state  is 
the  only  equilibrium  state  and  is  stable.  For  b  >  (3/2)77, 
the  isotropic  state  is  found  to  be  unstable.  The  first  nematic 
branch  emerging  at  b  =  (3/2)77  has  a  free  energy  lower 
than  that  of  the  isotropic  state,  and  thus  is  stable.  All  sub¬ 
sequent  nematic  branches  are  found  to  be  unstable  when 
perturbed  by  cos  29,  the  leading  Fourier  mode  in  the  first 
nematic  branch.  The  fundamental  difference  between  the 
Onsager  potential  and  the  Maier-Saupe  potential  is  that 
for  /  >  la  perturbation  of  the  Fourier  mode  cos  2  jO  in 
the  polymer  orientation  distribution  will  lead  to  a  corre¬ 
sponding  Fourier  mode  in  the  Onsager  potential  while  it 
has  no  effect  at  all  on  the  Maier-Saupe  potential.  In  other 
words,  for  /  >  1  the  Onsager  potential  responds  to  the 
Fourier  mode  cos  2 jd  in  the  polymer  orientation  distribu¬ 
tion  while  the  Maier-Saupe  potential  does  not.  The  /-th 
nematic  branch  can  be  viewed  as  nematic  state  excited  by 
a  perturbation  of  cos  2  jO  from  the  isotropic  state.  Although 
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all  nematic  branches  with  j  >  1  are  unstable  when  per¬ 
turbed  by  cos  26  and  consequently  is  physically  not  observ¬ 
able,  the  response  of  the  Onsager  interaction  potential  to  the 
Fourier  mode  cos  46  in  the  polymer  orientation  distribution 
is  still  significant.  If  an  external  field  of  the  form  of  cos  46 
is  applied,  it  will  induce  an  amplitude  of  cos  46  in  the  ori¬ 
entation  distribution  and  the  amplitude  will  be  enhanced  by 
the  Onsager  interaction  potential. 

Acknowledgments:  This  work  was  partially  supported 
by  the  Air  Force  Office  of  Scientific  Research  under 
grant  F1ATA06313G003  and  by  the  National  Science 
Foundation. 

References 

1.  D.  R.  J.  Chillingworth,  E.  Vicente  Alonso,  and  A.  A.  Wheeler, 
J.  Phys.  34,  1393  (2001). 

2.  F.  Cocchini,  C.  Aratari,  and  G.  Marruucci,  Macromolecules  23,  4446 

(1990). 

3.  P.  Constantin.  I.  Kevrekidis,  and  E.  S.  Titi,  Arch.  Rat.  Mech.  Anal. 
174,  365  (2004). 

4.  P.  Constantin,  I.  Kevrekidis,  and  E.  S.  Titi.  Discrete  and  Continuous 
Dynamical  Systems  11,  101  (2004). 

5.  P.  Constantin  and  J.  Vukadinovic,  Nonlinearity  18,  441  (2005). 

6.  M.  Doi  and  S.  F.  Edwards,  The  Theory  of  Polymer  Dynamics, 
Oxford  University  Press  (1986). 

7.  V.  Faraoni,  M.  Grosso,  S.  Crescitelli.  and  P.  L.  Maffetone,  .7.  Rheol. 
43,  829  (1999). 

8.  I.  Fatkullin  and  V.  Slastikov,  Nonlinearity  18,  2565  (2005). 

9.  M.  G.  Forest.  R.  Zhou,  and  Q.  Wang,  Phys.  Rev.  E  66.  031712 

(2002). 

10.  M.  G.  Forest,  Q.  Wang,  and  R.  Zhou,  Rheol.  Acta  44,  80  (2004). 

11.  M.  G.  Forest,  R.  Zhou,  and  Q.  Wang,  Rheol.  Acta  43,  17 

(2004). 

12.  M.  G.  Forest,  R.  Zhou,  and  Q.  Wang,  Phys.  Rev.  Lett.  93,  088301 

(2004). 

13.  M.  G.  Forest,  R.  Zhou,  and  Q.  Wang,  SIAM  MMS  3,  853  (2005). 

14.  M.  G.  Forest,  R.  Zhou,  and  Q.  Wang,  J.  Non-Newtonian  Fluid  Mech. 
116,  183  (2004). 

15.  M.  G.  Forest  and  Q.  Wang,  Rheol.  Acta  42,  20  (2003). 

16.  M.  G.  Forest,  Q.  Wang,  R.  Zhou,  and  E.  Choate,  J.  Non-Newtonian 
Fluid  Mech.  118,  17  (2004). 


17.  M.  G.  Forest,  Q.  Wang,  and  R.  Zhou,  J.  Rheol.  51,1  (2007). 

18.  D.  Grecov  and  A.  D.  Rey,  Rheologica  Ada  42,  590  (2003). 

19.  A.  Gopinath,  L.  Mahadevan,  and  R.  C.  Armstrong.  Physics  of  Fluids 
18,  028102  (2006). 

20.  S.  Z.  Hess,  Z.  Naturforsch.  A  31A,  1034  (1976). 

21.  N.  Kuzuu  and  M.  Doi,  Phys.  Soc.  Jpn.  52,  3486  (1983). 

22.  N.  Kuzuu  and  M.  Doi,  Phys.  Soc.  Jpn.  53,  1031  (1984). 

23.  J.  H.  Lee,  M.  G.  Forest,  and  R.  Zhou,  Discrete  and  Continuous 
Dynamical  Systems-Series  B  6,  339  (2006). 

24.  H.  Liu,  H.  Zhang,  and  P.  Zhang,  Comm.  Math.  Sci.  3,  201  (2005). 

25.  C.  Luo,  H.  Zhang,  and  P.  Zhang,  Nonlinearity  18,  379  (2005). 

26.  P.  L.  Maffettone  and  G.  Marrucci,  J.  Non-Newtonian  Fluid  Mech. 
38,  273  (1991). 

27.  P.  L.  Maffettone  and  S.  Crescitelli.  J.  Rheol.  38,  1559  (1994). 

28.  P  L.  Maffettone  and  S.  Crescitelli,  J.  Non-Newtonian  Fluid  Mech. 
59,  73  (1995). 

29.  P  L.  Maffettone,  M.  Grosso,  M.  C.  Friedenberg,  and  G.  G.  Fuller, 
Macromolecules  29,  8473  (1996). 

30.  G.  Marrucci  and  P.  L.  Maffettone,  Macromolecules  22,  4446  (1989). 

31.  G.  Marrucci  and  P.  L.  Maffettone,  J.  Rheol.  34,  1217  (1990). 

32.  T.  Maruyama,  G.  G.  Fuller,  M.  Grosso,  and  P  L.  Maffetone.  J.  Non- 
Newtonian  Fluid  Mech.  76,  233  (1998). 

33.  L.  Onsager,  Ann.  N.  Y.  Acad.  Sci.  51,  627  (1949). 

34.  A.  D.  Rey  and  M.  M.  Denn,  Annual  Review  of  Fluid  Mechanics 
34,  233  (2002). 

35.  G.  Rienacker  and  S.  Hess,  Physica  A  267,  294  (1999). 

36.  H.  See,  M.  Doi,  and  R.  G.  Larson,  .7.  Chem.  Phys.  92,  792  (1990). 

37.  A.  N.  Semenov,  Sov.  Phys.  JETP  66,  321  (1983). 

38.  G.  Sgalari,  G.  L.  Leal,  and  J.  J.  Feng,  J.  Non-Newtonian  Fluid 
Mech.  102,  361  (2002). 

39.  E.  Vicente  Alonso,  A.  A.  Wheeler,  and  T.  J.  Sluckin,  Proc.  R.  Soc. 
Fond.  A  459,  195  (2003). 

40.  H.  Wang  and  H.  Zhou,  Phys.  Lett.  A  372,  3423  (2008). 

41.  Q.  Wang,  S.  Sircar,  and  H.  Zhou,  Comm.  Math.  Sci.  3,  605  (2005). 

42.  K.  S.  Yim,  G.  G.  Fuller,  A.  Datko,  and  C.  Eisenbach,  Macro¬ 
molecules  34,  6972  (2001). 

43.  H.  Yu  and  P.  Zhang,  J.  Non-Newtonian  Fluid  Mech.  141,  1 16  (2007). 

44.  A.  Zarnescu,  Nonlinearity  19.  1619  (2006). 

45.  H.  Zhou,  H.  Wang,  M.  G.  Forest,  and  Q.  Wang,  Nonlinearity 
18,  2815  (2005). 

46.  H.  Zhou,  H.  Wang,  Q.  Wang,  and  M.  G.  Forest,  Nonlinearity  20,  277 

(2007). 

47.  H.  Zhou  and  H.  Wang,  Comm.  Math.  Sci.  5,  113  (2007). 

48.  H.  Zhou,  H.  Wang,  and  Q.  Wang,  Discrete  and  Continuous  Dynam¬ 
ical  Systems-Series  B  7,  907  (2007). 

49.  H.  Zhou  and  H.  Wang,  Physics  of  Fluids  19,  103107  (2007). 


Received:  10  June  2008.  Accepted:  30  December  2008. 


18 


J.  Comput.  Theor.  Nanosci.  7,  1-18,  2010 


